#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>

double psi(double alpha, double x){
	return exp(-alpha*x*x);
}

double EL(double alpha, double x){
	return alpha + (0.5-2*alpha*alpha)*x*x;
}

/*
 * vmc1d_aveEL:		变分蒙特卡洛方法，返回<EL>值
 * doubla alpha:	参数
 * int n:		number of random walker
 * double a, b:		x　轴定义域, 所有random walker被限制在[a,b]内
 * double h:		随机行走步长
 * int N:		每个 random walker的markov链长度
 */
double vmc1d_aveEL(double alpha, int n, 
		double a, double b, double h, int N){

	srand(time(NULL));
	int i,j,count;
	double x, d, y, z, p;
	y=0;
	for(i=0;i<n;i++){
		x = (double)rand()/RAND_MAX * (b-a) + a;//初始位置
		z=0;
		count=0;
		for(j=0;j<N;j++){
			if( j > 0.2 *N ){
				z += EL(alpha, x);
				count ++;
			}
			d = ( (rand()%2)*2 - 1 )*h; // -h or h
			if( x+d > b || x+d < a ) continue;
			p = pow( psi(alpha, x+d), 2) / pow( psi(alpha,x), 2);
			if(p>1) x += d;
			else{
				if( (double)rand()/RAND_MAX < p ) x+=d;
			}
		}
		z /= count;
		y += z;
	}
	return y/n;
}

int main(){

	ofstream fp("vmc1dho_aveEL.txt");
	for(double alpha=0.4; alpha <=0.6; alpha += 0.01){
		fp << alpha << "\t" << vmc1d_aveEL(alpha, 1000, -3, 3, 0.2, 1000) <<endl;
	}
	fp.close();
	return 0;
}
